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A recently proposed dimensional reduction approach for studying synchronization in the Kuramoto model is 
employed to build optimal network topologies to favor or to suppress synchronization. The approach is based 
in the introduction of a collective coordinate for the time evolution of the phase locked oscillators, in the spirit 
of the Ott-Antonsen ansatz. We show that the optimal synchronization of a Kuramoto network demands the 
maximization of the quadratic function a/Lai, where at stands for the vector of the natural frequencies of the 
oscillators, and L for the network Laplacian matrix. Many recently obtained numerical results can be re-obtained 
analytically and in a simpler way from our maximization condition. A computationally efficient hill climb 
rewiring algorithm is proposed to generate networks with optimal synchronization properties. Our approach can 
be easily adapted to the case of the Kuramoto models with both attractive and repulsive interactions, and again 
many recent numerical results can be rederived in a simpler and clearer analytical manner. 

PACS numbers: 89.75.Fb, 05.45.Xt, 89.75.Hc 


I. INTRODUCTION 


ramoto oscillators 0 is to use the order parameter defined as 


Synchronization 0] is present in a myriad of natural and 
synthetic systems, ranging from metabolic processes in popu¬ 
lations of yeast (2] to power grids 0, and its abundant pres¬ 
ence has stimulated in the last decades a very active area of 
research. The Kuramoto model BQ has been used as one 
of the most versatile tools to understand the different scenar¬ 
ios where a population of heterogeneous units can develop a 
common rhythm through mutual interaction, despite the in¬ 
trinsic tendency of each element to oscillate with its own nat¬ 
ural frequency when uncoupled from the network. In general, 
the units in question are located on the vertices of a complex 
undirected network 0, which is described by its adjacency 
matrix A, ;, with binary entries A (/ = 1 if there is an edge be¬ 
tween vertices i and j and A = 0 otherwise. The Kuramoto 
model is then defined by the nonlinear system of differential 
equations (9| 


dOj vn 

— = a>i + A 2^ A t j sin (0j - 6i), 
7=1 


( 1 ) 


where the phase $,(/) corresponds to the state of the ;th unit, 
which would oscillate with its natural frequency w, if uncou¬ 
pled from the network. The interaction between connected 
units is governed by the sine of their phase difference, and 
the strength of the coupling is determined by the parameter A. 
Typically, in order to mimic the inherent differences between 
the elements of real systems, the natural frequencies w,- are as¬ 
sumed to be random variables with a distribution g(a>), which 
we will consider here to be symmetric and unimodal. 

A convenient way to describe the global state of the Ku- 
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( 2 ) 


where r is assumed to be real, which corresponds to the 
centroid of the phases if they are considered as a swarm of 
points moving around the unit circle. For incoherent mo¬ 
tion, the phases are scattered on the circle homogeneously and 
r ~ N l/2 , while for a synchronized state they should move in 
a single lump and, consequently, r ~ . The Kuramoto model 
is known to exhibit a second order phase transition from the 
incoherent to the synchronized regime at a critical value A c of 
the coupling strength. For A > A, , the order parameter r is an 
increasing function of A. A very natural question to set forth 
here is: given N vertices with natural frequencies {cu,j^| and 
m edges, what is the best way to connect the vertices (form¬ 
ing a connected network) in order to optimize synchronization 
(maximize r)l In the last years, it has become clear mm 
that optimal (in the synchronization sense) networks of Ku¬ 
ramoto oscillators typically have a strong negative correlation 
between the natural frequencies of adjacent oscillators, and a 
positive correlation between the frequency magnitude |oj,| and 
the degree k, = Yl] \ A/, of the vertex where the i oscillator 
lies. This situation favors some heterogeneity in the degree 
distribution: vertices with large positive (negative) natural fre¬ 
quencies tend to have higher degrees and to be surrounded by 
vertices with negative (positive) natural frequencies. These 
properties are also observed in in the optimized networks 
for the Kuramoto model with inertia fl5l . although in this 
case the optimization has more severe consequences, as it also 
changes the nature of the phase transition: optimized networks 
typically possess a first order transition instead of the usual 
second order one. 

We notice also that many recent works have been devoted 
to the optimization of the synchronization for the case of iden¬ 
tical oscillators with diffusive coupling in a complex network. 
For such a case, the configuration where all the oscillators be¬ 
have identically is a solution of the equations of motion, al- 
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beit not necessarily a stable one. In order to guarantee the 
stability of this synchronized state Dana, one finds that 
the ratio R = w, v /tu 2 between the largest eigenvalue m N and 
the smallest non-zero eigenvalue m 2 of the Laplacian matrix 
L,j = djjkj - A;j must be smaller than a certain value, which, 
on the other hand, is uniquely determined by the internal dy¬ 
namics of the vertices (independent of the network topology, 
see also iT&HZOl ). Moreover, the smaller the ratio R, the bet¬ 
ter is the stability of the synchronized state. There has been a 
great amount of work in studying the properties of networks 
obtained from the minimization of the ratio R trough some 
kind of rewiring |211 , including in the case of multiplex net¬ 
works ED- We note, however, that these approaches are not 
suitable when the natural frequencies are random variables as 
in the cases considered here. 

In this work, we revisit the issue of optimization of syn¬ 
chronization in Kuramoto networks, but employing the di¬ 
mensional reduction approach recently proposed by Gottwald 
m , which explores some tools of the theory of solitary waves 
by introducing a collective coordinate for the time evolution of 
the phase locked oscillators, in the spirit of the Ott-Antonsen 
ansatz m ns. Thanks to this dimensional reduction, we 
are able to derive analytically a simple condition for optimiz¬ 
ing the network topology to favor synchronization in the Ku¬ 
ramoto model 0- We will show that the optimal synchro¬ 
nization of 0 for A > A c demands the minimization of the 
quadratic function 

N N N 

/'(0) = ^ 2 A ijUiWj - 'Yj k i^ = -u T Lu, (3) 

i= 1 7=1 i= 1 


where a> is the vector formed by the natural frequencies of 
the oscillators, and L the Laplacian matrix. Many previously 
known results can be readily obtained from the minimization 
of 0 or, in other words, from the maximization of lj t L(j). 
Moreover, our condition can be adapted into a hill climb al¬ 
gorithm to produce optimal networks in a very efficient way, 
since it is not necessary the integration of any differential 
equation or the computation of any matrix eigenvector and our 
condition involves only only plain matrix multiplications. We 
will also apply Gottwald’s dimensional reduction approach 
l23l to study the optimal networks for the Kuramoto model 
with attractive and repulsive interactions li26l 


dd, 

dt 


— OJj + 


YjAijAjSmiej-OiX 


(4) 


7=1 


where the coupling strengths Aj can now be either positive 
or negative, encoding not only the strength, but also the sign 
of the influence of oscillator j on its neighbors. Positive 
values of Aj promote in-phase relationships between neigh¬ 
bors, whereas negative values do anti-phase ones. For this 
kind of model, we also show that optimal synchronization de¬ 
mands the maximization of a quadratic function involving the 
Laplacian matrix associated to the weighted adjacency matrix 
Ch = A A , with A ij = djjAj. Our results in this case are also 
compatible with those ones recently discussed in Ii27ll29l . 

In the next section, we will derive our results for the Ku¬ 
ramoto model 0 , introduce our main algorithm, and discuss 


an explicit example. Section III is devoted to the extension of 
our approach to Kuramoto models with attractive and repul¬ 
sive interactions (0, while the last section is left to some final 
remarks. 


II. OPTIMAL SYNCHRONIZATION IN THE KURAMOTO 
MODEL 


Since the Kuramoto model 0 has rotational invariance, we 
can change, without loss of generality, to a new reference 
frame Q, —> 6, + Of such that the distribution g(a>) is cen¬ 
tered at co — 0 (which, in our case, implies (ai) = 0). Several 
well known results suggest that for the Kuramoto model 0, 
the time evolution of the phase locked oscillators may be well 
approximated by using the collective parametrization ||23l 

0,(t) = a(t)a>i. (5) 

By using such ansatz in the equations of motion 0 in the 
reference frame for which ( a >) = 0, and demanding that a 
certain error function is minimal, Gottwald end up with a one¬ 
dimensional differential equation for a l23l 

a = ( 6 ) 

cr z 

with cr 2 = and 


2 N N 

f(a) = — + ^ on ^ A tJ sin(a(wj - «,)). (7) 

*=1 J = 1 

With the ansatz 0, synchronized solutions correspond to 
stable fixed points a* of 0, as the difference of phases 
a(t ) (ojj - cu/j are constants if a = 0. Hence, the fixed points 
a* must satisfy the condition f(a*) = 0, and its stability is de¬ 
termined by the sign of Figure 0 shows f(a) for two 



FIG. 1. (Color online). The graphics show the function f(a) 0 for 
an Erdos-Renyi network with N = 1000 vertices, mean degree (k) = 
10 and natural frequencies drawn from the uniform distribution in 
the interval (—1,1), for two different values of the coupling strength 
A. 

different values of A, A —> co and A = 0.6, for an Erdos-Renyi 
network with N = 1000 vertices, mean degree (k) = 10 and 
natural frequencies drawn from the uniform distribution in the 
interval (-1,1). For the first case, we clearly see that there is 
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a stable fixed point a* = 0, as one would expect. The curve 
for A = 0.6, however, is the same one for A —» oo, but shifted 
upwards, displacing the stable fixed points towards right. 

This scenario allows us to devise a method to optimize the 
synchronization for finite values of A. We want to have r as 
large as possible for a given finet value of ,1 > A, . In order to 
increase the order parameter r, we must assure that the phase’s 
difference in the synchronized state are as small as possible, 
i.e, a * must be minimal, so that every oscillator is close to 
each other accordingly to ([5]). If we linearize the function f(ai) 
around a = 0 and demand that it crosses the abscissa as close 
as possible to the origin, we obtain the condition that /'(0) 
must be minimal, which, according to (|3]l implies that o/Loj 
must be maximal. Notice that, interestingly, this conditions 
does not depend on the value of A. It is clear also that optimal 
synchronization for the Kuramoto model demands a negative 
correlation of the natural frequencies of adjacent oscillators 
and a positive correlation between the degrees kj and the val¬ 
ues of \u>i\, which are precisely the results obtained previously 
by different and more intricate numerical approaches QSHED- 

A hill climb rewiring algorithm can be easily set up to find 
the maximum of lj t L u> and, hence, to produce optimal syn¬ 
chronization networks. The strategy is roughly the following. 
A randomly selected edge connecting two vertices is removed 
if it does not disconnect the network, and two randomly cho¬ 
sen disconnected vertices are then connected. After this step, 
the new value of u> T Lu> is evaluated. If the rewiring results in 
a higher value, one keeps the modification or, otherwise, one 
discharges the rewiring and returns the network to its previous 
configuration. This procedure is repeated until a> T Loj attains 
a minimum value, what of course will guarantee a minimal 
value for /'(0). In practice, our algorithm limits the maxi¬ 
mum number of iterations (up to 2 x 10 4 times in our simula¬ 
tions). These edge swaps preserve the average degree of the 
initial network since the number of edges is kept the same, 
but not the degree distribution. It is clear that this kind of hill 
climb algorithm works by performing a local search of the op¬ 
timal state by incrementally changing the stmcture of the net¬ 
work, and it is indeed the simplest algorithm to perform the 
maximization of u) T La>. Compared with more complex algo¬ 
rithms such as simulated annealing If30l . for instance, our hill 
climb algorithm has proved to be much faster and extremely 
reliable. However, there are some potential issues, as for in¬ 
stance, the algorithm might get stuck in a local minimal state 
far from the global minimum. Nevertheless, our numerical 
simulations show that this is very rare and we can indeed get 
networks with pronounced enhancement of the synchroniza¬ 
tion capacity from this simple algorithm, which is indeed the 
most used one in the literature of optimization of synchroniza¬ 
tion in complex networks (see, for instance, 110-121). 

As an explicit example, we apply this algorithm for a net¬ 
work built with the Barabasi-Albert method Oil (which cre¬ 
ates scale free networks with degree distribution p(k) oc k 3 ) 
with N = 10 3 vertices and mean degree (k) = 6. The natu¬ 
ral frequencies were drawn from the unit normal distribution. 
The synchronization diagram for this network is shown as red 
squares in Fig. [2] where the usual smooth monotonically in¬ 
creasing behavior for r over all the range of coupling strength 



FIG. 2. (Color online). Synchronization diagrams for the initial and 
optimized networks for the Kuramoto model 0 - The diagram for 
the initial network, a Barabasi-Albert network with N = 10 3 vertices 
and ( k ) = 6, is the curve with the red squares. The natural frequen¬ 
cies were drawn from the unit normal distribution. The optimized 
network, shown as the curve with blue circles, was obtained from the 
maximization of lj t Lcj by using the hill climb algorithm described 
in the text. The inset shows the region around the phase transition for 
the rewired network, confirming that is of the second order type. 


A is observed. The synchronization diagram of the optimized 
network that results from the algorithm is shown in Figure [2] 
as blue circles. It is evident the considerable enhancement 
of the network synchronization capacity over a large range of 
coupling strengths. We stress that the phase transition is now 
much more abrupt, but still of second order, since the zoom 
around the critical region shown in the inset has a continuous 
behavior, and no hysteresis loop seems to be present. Figure 
[3] shows the positive correlation between natural frequencies 
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FIG. 3. (Color online). The graphics shows in the left panel the posi¬ 
tive correlation between natural frequencies |o»| of the oscillators and 
degrees k of the vertices where they lie. The negative correlation be¬ 
tween natural frequencies of adjacent vertices, measured here by the 
average natural frequency (w,)/v = Ayu^/k,- of the neighbours 
Nj of a vertex i as a function of its natural frequency at h is depicted 
in the right panel. 

\u)\ and degrees k for the optimized networks (left panel), as 
well as the negative correlation between the natural frequen¬ 
cies of adjacent oscillators (right panel). We can give yet an¬ 
other measure of this strong negative correlation noticing that 
the fraction of links connecting oscillators with positive and 
negative frequencies jumps from 0.5 to 0.87 in this particular 
run of the optimization procedure. 

The condition of maximal u> t Lcj allows us also to attack 
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another related problem considered recently in the literature. 
The Laplacian matrix L of an undirected network is symmet¬ 
ric, so that all its eigenvalues are real. Furthermore it can be 
demonstrated (8) that for a connected network it has only one 
eigenvalue tu\ — 0 and all the other N - \ are positive. Sup¬ 
pose that the topology of the network is fixed and we want 
to obtain the set of natural frequencies to be allocated to each 
vertex such that it optimizes the synchronization properties, 
subjected to the constraint cr 2 = w 2 = N (in fact, any 
positive constant). By writing w as a linear combination of 
the eigenvalues V; of L, u> — Yi?=i Q i v i, where of = N to 
assure the validity of constraint, we have that in order to max¬ 
imize u) T Llj one must set co proportional to the eigenvector of 
largest eigenvalue of L (a i = ... = a^-\ = 0 and a n = V N), a 
result recently found in Ifl3l by using more intricate methods. 


In the same manner as for the Kuramoto model, to obtain the 
stable synchronized solutions we must seek values of a* such 
that h(a*) = 0 with h'(a*) < 0. Two examples of the func¬ 
tion li(a) are shown in Figure [4] for an Erdos-Renyi network 



III. THE KURAMOTO MODEL WITH ATTRACTIVE AND 
REPULSIVE INTERACTIONS 


For the Kuramoto model with attractive and repulsive inter¬ 
actions 0. we will take advantage of the rather unexpected 
result that for both kinds of oscillator, those with positive and 
negative couplings, the collective parametrisation of the phase 
locked oscillators can be taken to be the same, see, for in¬ 
stance, Figure 3 of ll26l . Hence, we can employ the same 
ansatz 0 also in this case. We will consider the situation 
with two kinds of units: the attractive and repulsive ones, for 
which the coupling strength are, respectively, A 1 and A . Let 
p be the fraction of attractive oscillators, and, consequently, 
1 — p the fraction of repulsive ones. As the coupling strengths 
are properties of each oscillator, it is possible that an oscillator 
with positive coupling strength is connected to an oscillator 
with negative strength, causing in this way a kind of frustra¬ 
tion effect in the network. Without loss of generality, one can 
rescale the time variable in order to have T + = 1. An issue 
that arises here is that the frequency of rotation of the locked 
oscillators is not equal to the average value of the natural fre¬ 
quencies ( u >), as it happens for the Kuramoto model ([TJ. By 
changing to a new reference frame 0, —> 0, + Of, one has from 

0 


dO, 

dt 


N 

= ojj - Q + y AjjAj sin(0y - 0,). (8) 

j= i 


Multiplying both sides by A, and then summing with respect 
to the index i, the term 2^, Z yii ^ijAiAj sin(0y - 0,) cancel 
out, and we have that a synchronized state must have 


Zh AjCJi 

Zf=i ^ ' 


(9) 


Repeating here the same steps of Gottwald lf23l . we will have 
that a(t ) must obey, in this case, the ordinary differential equa¬ 
tion a = h(a ), where 


1 N N 

h(a ) = 1 4-- ^ 0J i ^ AijAj sm(a(cjj - cu,)). (10) 

<T i= i j= i 


FIG. 4. (Color online). The graphics show the function h(a) ( | 10[ > for 
an Erdos-Renyi network with N = 10 3 vertices and natural frequen¬ 
cies drawn from the unit normal distribution, for two different values 
of the fraction p of vertices with positive coupling strength. For both 
curves, A~ = -0.5. 


with N = 10 3 vertices, mean degree (k) = 10 and natural fre¬ 
quencies drawn from a unit normal distribution. The curves 
correspond to two different values of the fraction p of oscilla¬ 
tors with positive coupling strengths (and ,1 = -0.5). From 
the overall shape of the function h(a), the condition for the 
optimization, either to favor or to suppress synchronization, 
depends on the value of the derivative of li(a) at a - 0, 

N N N N 

O-V(O) = Y^^AijAjCOiUjj - Y j Yj A 0 A j oj1 i = 

i=l 7=1 '=1 7=1 

(ID 

where £. is the Laplacian matrix associated to the weighted 
adjacency matrix Aft = AA, 


-£(/ — ^ij ^ik dftij, 


( 12 ) 


k=l 


with A ij = 6ijAj. Clearly, when T, = A for all oscillators, we 
recover the previous optimization condition. In contrast with 
the previous case, the Laplacian matrix X is not symmetrical 
here, rendering the spectral analysis much more complicate in 
this case. 

However, the overall scenario is very similar to that one 
found for the standard Kuramoto model, and we will proceed 
in the same manner to optimize the network structure either to 
favor or to suppress synchronization. Of course, for the first 
case, we must maximize w r -Cw, whereas for the second case 
we must minimize it. It is interesting that we also have here 
a rather clear interpretation of the characteristics of the opti¬ 
mized network in terms of correlations between microscopic 


properties. Inspecting equation (11 1 , we see that if we are op¬ 


timizing the network to favor synchronization, there must be 
a negative correlation between the natural frequency co, of the 
oscillator placed at vertex i and the average value of the prod¬ 
uct ojjAj over its neighbors, as well as a positive correlation 
between |<y,j and the average value of coupling strengths of 
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its neighbors. If we are optimizing the network to suppress 
synchronization, the correlations must be reversed. The same 
hill climb rewiring algorithm can be implemented to minimize 
or maximize condition cj t £.u> and produce networks with the 
desired properties. 

We apply our algorithm to build the optimal network 
topologies, both to favor and to suppress synchronization, us¬ 
ing as seed an Erdos-Renyi network with 10 3 vertices and 
mean degree ( a >) = 5. The natural frequencies were drawn 
from the unit normal distribution and a fraction p = 0.8 of 
the vertices have positive strengths d + = 1, while 0.2 of the 
vertices have A~ = -0.5. The natural frequencies and cou¬ 
pling strengths for each vertex were randomly assigned and 
were kept fixed during the optimization procedure, the algo¬ 
rithm only performs the rewiring of the network connections. 
The results are shown in Figure [5] Panel (a) shows the time 



u>j of the oscillator at vertex i and the mean value 


(wT)jv = 


Z /_ i A-ijOtjAj 

Z.%i A u*j 


(13) 


over its set of neighbors JV,, and the positive correlation be¬ 
tween |w;| and the mean value of the coupling (A) of its neigh¬ 
bors JV,-. We have found also a positive correlation between 
the magnitude of the natural frequencies and degrees for the 
optimized network to favor synchronization, see Figure [6] 
whereas no correlation (either positive or negative) is seem 



FIG. 6. (Color online). The graphics show the positive correlation 
between the magnitude of the natural frequencies and degrees for the 
optimized network to favor synchronization, panel (a), whereas there 
is no correlation for the network optimized to suppress synchroniza¬ 
tion, panel (b). Both results correspond to the networks depicted in 
figure [5] The insets show, for both cases, the distribution of the cou¬ 
plings as function of the degrees. 


FIG. 5. (Color online). Panel (a) shows the time evolution of r(t) 
for the original network (the red dashed line), for the optimized net¬ 
work in order to favor synchronization (the blue continuous line) and 
for the optimized network in order to suppress synchronization (the 
dash-dotted green line). Panel (b) shows the behavior of function 
h(a) for the three cases depicted in panel (a). Panels (c) and (d) show, 
respectively, the negative correlation between the natural frequency 
w, of the oscillator at vertex i and the mean value of the product (a>A) 
over its set of neighbors )V, and the positive correlation between cjj 
and the mean value of the coupling (d) of (V, for the optimized net¬ 
work to favor synchronization shown in panel (a). 

evolution of the order parameter r(t) for the original network 
(the red dashed line), for the optimized network in order to 
favor synchronization (the blue continuous line) and for the 
optimized network in order to suppress synchronization (the 
dash-dotted green line). The overall behavior of the three 
cases are independent of the initial conditions. The method 
works very well for both cases, the synchronization proper¬ 
ties are clearly enhanced for the optimal synchronization net¬ 
work, while the network optimized to suppress synchroniza¬ 
tion shows only a noise-like behavior of the order parameter r. 
The panel (b) shows the functions h(a) for the three cases de¬ 
picted in panel (a). We can see that for the optimized network 
to suppress synchronization, the derivative h'( 0) reversed its 
sign, destroying the stable synchronized solution. Panels (c) 
and (d) show, for the network optimized to favor synchroniza¬ 
tion, the presence of the correlations discussed earlier, respec¬ 
tively, the negative correlation between the natural frequency 


for the optimized network to suppress synchronization. An¬ 
other interesting result is that for the optimal synchronization 
case, the positive and negative coupling strengths are equally 
distributed in a large range of degrees, while for the case of 
optimized network to suppress synchronization, most of the 
oscillators with negative coupling are located at large degree 
vertices (see the insets of Figure|6|. Notice that the optimized 
network to favor synchronization does not show r ~ 1, it in¬ 
stead saturates at the lower value r ~ 0.9 due to the already 
commented frustrating effects associated with the combina¬ 
tion of attractive and negative units. 

IV. CONCLUSIONS 

We have shown that the dimensional reduction approach re¬ 
cently proposed by Gottwald li23l for the Kuramoto model can 
be adapted to obtain a simple analytical condition to optimize 
the network topology to favor or to suppress synchronization. 
Our condition allowed us to rederive analytically some recent 
results mm obtained numerically. Our approach was also 
extended to the Kuramoto model with attractive and repulsive 
interactions, leading also to simple analytical condition to op¬ 
timize the networks, complementing in this way some recent 
works I !27]j29ll which considered the optimization of synchro¬ 
nization in these models numerically. In both cases, the opti¬ 
mal synchronization condition corresponds to the maximiza¬ 
tion of the quadratic funcion lj t L<jj, where where oj stands for 
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the vector of the natural frequencies of the oscillators, and L 
for the pertinent network Laplacian matrix. The optimization 
condition involves only the microscopic properties of the net¬ 
work, enlightening in this way many correlations observed nu¬ 
merically for optimal synchronization networks. We could in¬ 
troduce a hill climb rewiring algorithm to produce optimized 
networks in a very efficient way, since it is not necessary the 
integration of any differential equation or the computation of 
any matrix eigenvector, only plain matrix multiplications are 
used in each step. As common for this kind of algorithm, one 
cannot assure in principle that the procedure will stop effec¬ 
tively in a global minimum. However, since the algorithm is 
computationally simple, on can run the procedure for several 
initial networks in order to seek for a global minimum. In 
all tests we have performed, our algorithm returned, with lit¬ 
tle computational effort, networks with greatly enhanced syn¬ 
chronization properties. An important open problem is the 


optimization of other kinds of oscillators as, for example, the 
Winfree model l32l . It would be interesting to investigate if 
the Gottwald approach 1(231 . or some variant of it, could be 
applied to these cases, and if an analysis analogous to the 
presented here could be indeed used to optimize the network 
topology to favor synchronization or other kinds of collective 
behaviors that might emerge in these more complicated mod¬ 
els. 
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